Практическое задание к уроку 1 (2 неделя).

Линейная регрессия: переобучение и регуляризация

В этом задании мы на примерах увидим, как переобучаются линейные модели, разберем, почему так происходит, и выясним, как диагностировать и контролировать переобучение.

Во всех ячейках, где написан комментарий с инструкциями, нужно написать код, выполняющий эти инструкции. Остальные ячейки с кодом (без комментариев) нужно просто выполнить. Кроме того, в задании требуется отвечать на вопросы; ответы нужно вписывать после выделенного слова "Ответ:".

Напоминаем, что посмотреть справку любого метода или функции (узнать, какие у нее аргументы и что она делает) можно с помощью комбинации Shift+Tab. Нажатие Tab после имени объекта и точки позволяет посмотреть, какие методы и переменные есть у этого объекта.


In [3]:
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

Мы будем работать с датасетом "bikes_rent.csv", в котором по дням записаны календарная информация и погодные условия, характеризующие автоматизированные пункты проката велосипедов, а также число прокатов в этот день. Последнее мы будем предсказывать; таким образом, мы будем решать задачу регрессии.

Знакомство с данными

Загрузите датасет с помощью функции pandas.read_csv в переменную df. Выведите первые 5 строчек, чтобы убедиться в корректном считывании данных:


In [4]:
# (0 баллов)
# Считайте данные и выведите первые 5 строк
df = pd.read_csv('bikes_rent.csv')
df.head(5)


Out[4]:
season yr mnth holiday weekday workingday weathersit temp atemp hum windspeed(mph) windspeed(ms) cnt
0 1 0 1 0 6 0 2 14.110847 18.18125 80.5833 10.749882 4.805490 985
1 1 0 1 0 0 0 2 14.902598 17.68695 69.6087 16.652113 7.443949 801
2 1 0 1 0 1 1 1 8.050924 9.47025 43.7273 16.636703 7.437060 1349
3 1 0 1 0 2 1 1 8.200000 10.60610 59.0435 10.739832 4.800998 1562
4 1 0 1 0 3 1 1 9.305237 11.46350 43.6957 12.522300 5.597810 1600

Для каждого дня проката известны следующие признаки (как они были указаны в источнике данных):

  • season: 1 - весна, 2 - лето, 3 - осень, 4 - зима
  • yr: 0 - 2011, 1 - 2012
  • mnth: от 1 до 12
  • holiday: 0 - нет праздника, 1 - есть праздник
  • weekday: от 0 до 6
  • workingday: 0 - нерабочий день, 1 - рабочий день
  • weathersit: оценка благоприятности погоды от 1 (чистый, ясный день) до 4 (ливень, туман)
  • temp: температура в Цельсиях
  • atemp: температура по ощущениям в Цельсиях
  • hum: влажность
  • windspeed(mph): скорость ветра в милях в час
  • windspeed(ms): скорость ветра в метрах в секунду
  • cnt: количество арендованных велосипедов (это целевой признак, его мы будем предсказывать)

Итак, у нас есть вещественные, бинарные и номинальные (порядковые) признаки, и со всеми из них можно работать как с вещественными. С номинальныеми признаками тоже можно работать как с вещественными, потому что на них задан порядок. Давайте посмотрим на графиках, как целевой признак зависит от остальных


In [3]:
fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(15, 10))
for idx, feature in enumerate(df.columns[:-1]):
    df.plot(feature, "cnt", subplots=True, kind="scatter", ax=axes[idx / 4, idx % 4])


Блок 1. Ответьте на вопросы (каждый 0.5 балла):

  1. Каков характер зависимости числа прокатов от месяца?
    • ответ: Похожа на - квадратичную, то есть к лету достигает максимума, постепенно увеличивая значения, потом после лета постепенно падает
  2. Укажите один или два признака, от которых число прокатов скорее всего зависит линейно
    • ответ: от temp и от atemp

Давайте более строго оценим уровень линейной зависимости между признаками и целевой переменной. Хорошей мерой линейной зависимости между двумя векторами является корреляция Пирсона. В pandas ее можно посчитать с помощью двух методов датафрейма: corr и corrwith. Метод df.corr вычисляет матрицу корреляций всех признаков из датафрейма. Методу df.corrwith нужно подать еще один датафрейм в качестве аргумента, и тогда он посчитает попарные корреляции между признаками из df и этого датафрейма.


In [5]:
# Код 1.1 (0.5 балла)
# Посчитайте корреляции всех признаков, кроме последнего, с последним с помощью метода corrwith:
df[df.columns[:-1]].corrwith(df["cnt"])


Out[5]:
season            0.406100
yr                0.566710
mnth              0.279977
holiday          -0.068348
weekday           0.067443
workingday        0.061156
weathersit       -0.297391
temp              0.627494
atemp             0.631066
hum              -0.100659
windspeed(mph)   -0.234545
windspeed(ms)    -0.234545
dtype: float64

В выборке есть признаки, коррелирующие с целевым, а значит, задачу можно решать линейными методами.

По графикам видно, что некоторые признаки похожи друг на друга. Поэтому давайте также посчитаем корреляции между вещественными признаками.


In [5]:
# Код 1.2 (0.5 балла)
# Посчитайте попарные корреляции между признаками temp, atemp, hum, windspeed(mph), windspeed(ms) и cnt
# с помощью метода corr:
df[['temp','atemp', 'hum','windspeed(mph)','windspeed(ms)','cnt']].corr()


Out[5]:
temp atemp hum windspeed(mph) windspeed(ms) cnt
temp 1.000000 0.991702 0.126963 -0.157944 -0.157944 0.627494
atemp 0.991702 1.000000 0.139988 -0.183643 -0.183643 0.631066
hum 0.126963 0.139988 1.000000 -0.248489 -0.248489 -0.100659
windspeed(mph) -0.157944 -0.183643 -0.248489 1.000000 1.000000 -0.234545
windspeed(ms) -0.157944 -0.183643 -0.248489 1.000000 1.000000 -0.234545
cnt 0.627494 0.631066 -0.100659 -0.234545 -0.234545 1.000000

На диагоналях, как и полагается, стоят единицы. Однако в матрице имеются еще две пары сильно коррелирующих столбцов: temp и atemp (коррелируют по своей природе) и два windspeed (потому что это просто перевод одних единиц в другие). Далее мы увидим, что этот факт негативно сказывается на обучении линейной модели.

Напоследок посмотрим средние признаков (метод mean), чтобы оценить масштаб признаков и доли 1 у бинарных признаков.


In [6]:
# Код 1.3 (0.5 балла)
# Выведите средние признаков
print df.mean()


season               2.496580
yr                   0.500684
mnth                 6.519836
holiday              0.028728
weekday              2.997264
workingday           0.683995
weathersit           1.395349
temp                20.310776
atemp               23.717699
hum                 62.789406
windspeed(mph)      12.762576
windspeed(ms)        5.705220
cnt               4504.348837
dtype: float64

Признаки имеют разный масштаб, значит для дальнейшей работы нам лучше нормировать матрицу объекты-признаки.

Проблема первая: коллинеарные признаки

Итак, в наших данных один признак дублирует другой, и есть еще два очень похожих. Конечно, мы могли бы сразу удалить дубликаты, но давайте посмотрим, как бы происходило обучение модели, если бы мы не заметили эту проблему.

Для начала проведем масштабирование, или стандартизацию признаков: из каждого признака вычтем его среднее и поделим на стандартное отклонение. Это можно сделать с помощью метода scale.

Кроме того, нужно перемешать выборку, это потребуется для кросс-валидации.


In [9]:
from sklearn.preprocessing import scale
from sklearn.utils import shuffle

In [10]:
df_shuffled = shuffle(df, random_state=123)
X = scale(df_shuffled[df_shuffled.columns[:-1]])
y = df_shuffled["cnt"]

Давайте обучим линейную регрессию на наших данных и посмотрим на веса признаков.


In [6]:
from sklearn.linear_model import LinearRegression

In [11]:
# Код 2.1 (1 балл)
# Создайте объект линейного регрессора, обучите его на всех данных и выведите веса модели 
# (веса хранятся в переменной coef_ класса регрессора).
# Можно выводить пары (название признака, вес), воспользовавшись функцией zip, встроенной в язык python
# Названия признаков хранятся в переменной df.columns
alg_lnr = LinearRegression()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)


[('season', 570.87160362866916), ('yr', 1021.9663731167183), ('mnth', -141.30760553455309), ('holiday', -86.758235317525887), ('weekday', 137.22356194502672), ('workingday', 56.39073082795796), ('weathersit', -330.22867178815972), ('temp', 367.46506778728309), ('atemp', 585.56759062259221), ('hum', -145.61127178742132), ('windspeed(mph)', 12458844098749.342), ('windspeed(ms)', -12458844098947.797)]

Мы видим, что веса при линейно-зависимых признаках по модулю значительно больше, чем при других признаках.

Чтобы понять, почему так произошло, вспомним аналитическую формулу, по которой вычисляются веса линейной модели в методе наименьших квадратов:

$w = (X^TX)^{-1} X^T y$.

Если в X есть коллинеарные (линейно-зависимые) столбцы, матрица $X^TX$ становится вырожденной, и формула перестает быть корректной. Чем более зависимы признаки, тем меньше определитель этой матрицы и тем хуже аппроксимация $Xw \approx y$. Такая ситуацию называют проблемой мультиколлинеарности, вы обсуждали ее на лекции.

С парой temp-atemp чуть менее коррелирующих переменных такого не произошло, однако на практике всегда стоит внимательно следить за коэффициентами при похожих признаках.

Решение проблемы мультиколлинеарности состоит в регуляризации линейной модели. К оптимизируемому функционалу прибавляют L1 или L2 норму весов, умноженную на коэффициент регуляризации $\alpha$. В первом случае метод называется Lasso, а во втором --- Ridge. Подробнее об этом также рассказано в лекции.

Обучите регрессоры Ridge и Lasso с параметрами по умолчанию и убедитесь, что проблема с весами решилась.


In [12]:
from sklearn.linear_model import Lasso, Ridge

In [12]:
# Код 2.2 (0.5 балла)
# Обучите линейную модель с L1-регуляризацией
alg_lnr = Lasso()
alg_lnr.fit(X, y)

print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)


[('season', 560.24161603088692), ('yr', 1019.4634940657198), ('mnth', -128.73062703678767), ('holiday', -86.152781333711047), ('weekday', 137.34789390496346), ('workingday', 55.21237064135687), ('weathersit', -332.36985696234876), ('temp', 376.36323620969893), ('atemp', 576.53079350455039), ('hum', -144.12915500348589), ('windspeed(mph)', -197.1396894024866), ('windspeed(ms)', -2.8049200120791778e-08)]

In [13]:
# Код 2.3 (0.5 балла)
# Обучите линейную модель с L2-регуляризацией
alg_lnr = Ridge()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)


[('season', 563.064572252007), ('yr', 1018.948378787532), ('mnth', -131.87332028247087), ('holiday', -86.74609799709188), ('weekday', 138.00511117871915), ('workingday', 55.903110375064927), ('weathersit', -332.34978849907424), ('temp', 386.45788919200811), ('atemp', 566.34704705996865), ('hum', -145.07132729867027), ('windspeed(mph)', -99.259441081952488), ('windspeed(ms)', -99.259441154183193)]

Блок 2. Поясните, каким образом введение регуляризации решает проблему с весами и мультиколлинеарностью.

Ваш ответ (1 балл): регуляризация изменяет функцию ошибку, таким образом модель "наказывается" за слишком большие коэф.. Коэф. перед скоростью там получился просто коллосальный на много порядков больше, чем перед остальными, регуляризация же решила эту проблему

Проблема вторая: неинформативные признаки

В отличие от L2-регуляризации, L1 обнуляет веса при некоторых признаках. Объяснение данному факту дается в одной из лекций курса.

Давайте пронаблюдаем, как меняются веса при увеличении коэффициента регуляризации $\alpha$ (в лекции коэффициент при регуляризаторе мог быть обозначен другой буквой).


In [14]:
# Код 3.1 (1 балл)
alphas = np.arange(1, 500, 50)
coefs_lasso = np.zeros((alphas.shape[0], X.shape[1])) # матрица весов размера (число регрессоров) x (число признаков)
coefs_ridge = np.zeros((alphas.shape[0], X.shape[1]))
# Для каждого значения коэффициента из alphas обучите регрессор Lasso
# и запишите веса в соответствующую строку матрицы coefs_lasso (вспомните встроенную в python функцию enumerate),
# а затем обучите Ridge и запишите веса в coefs_ridge.
print coefs_lasso.shape
num = 0
for alpha  in alphas:
    modelL1 = Lasso(alpha = alpha).fit(X, y)
    modelL2 = Ridge(alpha = alpha).fit(X, y)

    coefs_lasso[num] = modelL1.coef_
    coefs_ridge[num] = modelL2.coef_
    
    num = num + 1


(10, 12)

Визуализируем динамику весов при увеличении параметра регуляризации:


In [15]:
plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_lasso.T, df.columns):
    plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Lasso")

plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_ridge.T, df.columns):
    plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Ridge")

plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_lasso.T, df.columns):
    if feature == 'windspeed(mph)' or feature == 'windspeed(ms)' or feature == 'yr':
        plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Lasso - temp")


Out[15]:
<matplotlib.text.Text at 0xca3ee70>

Ответы на следующие вопросы можно давать, глядя на графики или выводя коэффициенты на печать.

Блок 3. Ответьте на вопросы (каждый 0.25 балла):

  1. Какой регуляризатор (Ridge или Lasso) агрессивнее уменьшает веса при одном и том же alpha?
    • Ответ:Lasso
  2. Что произойдет с весами Lasso, если alpha сделать очень большим? Поясните, почему так происходит.
    • Ответ: многие из них станут = 0, такое происходит в силу особенности функционала ошибки, при котором на ошибку влияет модуль сумму коэф.
  3. Можно ли утверждать, что Lasso исключает один из признаков windspeed при любом значении alpha > 0? А Ridge? Ситается, что регуляризатор исключает признак, если коэффициент при нем < 1e-3.
    • Ответ: Да можно, это видно на 3ем графике выше.
  4. Какой из регуляризаторов подойдет для отбора неинформативных признаков?
    • Ответ: L1, это можно определить по 0-ым коэф. перед ними

Далее будем работать с Lasso.

Итак, мы видим, что при изменении alpha модель по-разному подбирает коэффициенты признаков. Нам нужно выбрать наилучшее alpha.

Для этого, во-первых, нам нужна метрика качества. Будем использовать в качестве метрики сам оптимизируемый функционал метода наименьших квадратов, то есть Mean Square Error.

Во-вторых, нужно понять, на каких данных эту метрику считать. Нельзя выбирать alpha по значению MSE на обучающей выборке, потому что тогда мы не сможем оценить, как модель будет делать предсказания на новых для нее данных. Если мы выберем одно разбиение выборки на обучающую и тестовую (это называется holdout), то настроимся на конкретные "новые" данные, и вновь можем переобучиться. Поэтому будем делать несколько разбиений выборки, на каждом пробовать разные значения alpha, а затем усреднять MSE. Удобнее всего делать такие разбиения кросс-валидацией, то есть разделить выборку на K частей, или блоков, и каждый раз брать одну из них как тестовую, а из оставшихся блоков составлять обучающую выборку.

Делать кросс-валидацию для регрессии в sklearn совсем просто: для этого есть специальный регрессор, LassoCV, который берет на вход список из alpha и для каждого из них вычисляет MSE на кросс-валидации. После обучения (если оставить параметр cv=3 по умолчанию) регрессор будет содержать переменную mse_path_, матрицу размера len(alpha) x k, k = 3 (число блоков в кросс-валидации), содержащую значения MSE на тесте для соответствующих запусков. Кроме того, в переменной alpha_ будет храниться выбранное значение параметра регуляризации, а в coef_, традиционно, обученные веса, соответствующие этому alpha_.

Обратите внимание, что регрессор может менять порядок, в котором он проходит по alphas; для сопоставления с матрицей MSE лучше использовать переменную регрессора alphas_.


In [14]:
from sklearn.linear_model import LassoCV

In [18]:
# Код 3.2 (1 балл)
# Обучите регрессор LassoCV на всех параметрах регуляризации из alpha
# Постройте график _усредненного_ по строкам MSE в зависимости от alpha. 
# Выведите выбранное alpha, а также пары "признак-коэффициент" для обученного вектора коэффициентов
alphas = np.arange(1, 100, 5)

clf = LassoCV(alphas = alphas)
model = clf.fit(X, y)
print model.mse_path_.shape
plt.plot(model.alphas_, model.mse_path_.mean(axis =1))
plt.xlabel('alpha')
plt.ylabel('mean MSE')
plt.show()


print 'выбранное alpha = ', model.alpha_,'\n', zip(df.columns[0:(df.shape[1]-1)], clf.coef_)


(20, 3)
[ 851097.97556673  844888.21156583  839042.30891536  833559.17934129
  828438.72799191  823680.87676286  819285.55495756  815252.69915064
  811582.25286145  808274.1660866   805329.65107759  802745.81156384
  800525.55750255  798642.91265887  796928.77526934  795177.42985684
  793826.01832989  792028.11513622  790397.72811014  825435.81203957]
[ 851097.97556673  844888.21156583  839042.30891536  833559.17934129
  828438.72799191  823680.87676286  819285.55495756  815252.69915064
  811582.25286145  808274.1660866   805329.65107759  802745.81156384
  800525.55750255  798642.91265887  796928.77526934  795177.42985684
  793826.01832989  792028.11513622  790397.72811014  825435.81203957]
выбранное alpha =  6 
[('season', 532.01898284135382), ('yr', 1015.0602226430597), ('mnth', -100.03952614356666), ('holiday', -83.293959875299393), ('weekday', 132.50446549095801), ('workingday', 51.557085614073799), ('weathersit', -330.55985673998117), ('temp', 370.67985503003416), ('atemp', 581.39693106549805), ('hum', -140.00740550068912), ('windspeed(mph)', -191.77140847135016), ('windspeed(ms)', -2.6855109161225938e-08)]

Итак, мы выбрали некоторый параметр регуляризации. Давайте посмотрим, какие бы мы выбирали alpha, если бы делили выборку только один раз на обучающую и тестовую, то есть рассмотрим траектории MSE, соответствующие отдельным блокам выборки.


In [19]:
# Код 3.3 (1 балл)
# Выведите значения alpha, соответствующие минимумам MSE на каждом разбиении (то есть по столбцам).
# На трех отдельных графиках визуализируйте столбцы .mse_path_

print model.alphas_[np.argmin(model.mse_path_[:,0])]
plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,0], 'r')
plt.xlabel('alpha')
plt.ylabel('1st MSE')
plt.legend()
plt.axhline(y= np.amin(model.mse_path_[:,0],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,0])]))

print model.alphas_[np.argmin(model.mse_path_[:,1])]

plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,1], color = 'y')
plt.xlabel('alpha')
plt.ylabel('2nd MSE')
plt.axhline(y= np.amin(model.mse_path_[:,1],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,1])]))

print model.alphas_[np.argmin(model.mse_path_[:,2])]

plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,2], color = 'g')
plt.xlabel('alpha')
plt.ylabel('3rd MSE')
plt.axhline(y= np.amin(model.mse_path_[:,2],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,2])]))


41
6
1
C:\Users\SBT-Ashrapov-IR\AppData\Local\Continuum\Anaconda2\lib\site-packages\matplotlib\axes\_axes.py:519: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "
Out[19]:
<matplotlib.text.Text at 0x9f8b250>

На каждом разбиении оптимальное значение alpha свое, и ему соответствует большое MSE на других разбиениях. Получается, что мы настраиваемся на конкретные обучающие и контрольные выборки. При выборе alpha на кросс-валидации мы выбираем нечто "среднее", что будет давать приемлемое значение метрики на разных разбиениях выборки.

Наконец, как принято в анализе данных, давайте проинтерпретируем результат.

Блок 4. Ответьте на вопросы (каждый 0.5 балла):

  1. В последней обученной модели выберите 4 признака с наибольшими (положительными) коэфициентами (и выпишите их), посмотрите на визуализации зависимостей cnt от этих признаков, которые мы рисовали в блоке "Знакомство с данными". Видна ли возрастающая линейная зависимость cnt от этих признаков по графикам? Логично ли утверждать (из здравого смысла), что чем больше значение этих признаков, тем больше людей захотят взять велосипеды?
    • Ответ: yr, season, atemp, temp. Чем больше значение признака, тем больше cnt, но почему так - не так очевидно. Рост температуры и проката связан с тем, что просто зимой не покатаешься, летом же это уже не так важно. Осенью меньше чем летом прокат. А второй год больше значение в силу роста самого бизнеса. Не все так однозначно.
  2. Выберите 3 признака с наибольшими по модулю отрицательными коэффициентами (и выпишите их), посмотрите на соответствующие визуализации. Видна ли убывающая линейная зависимость? Логично ли утверждать, что чем больше величина этих признаков, тем меньше людей захотят взять велосипеды?
    • Ответ: weathersit, hum, windspeed. Да видна - очевидно при дожде и ветре кататься не хочется, как и при плохом прогнозе на погоду.
  3. Выпишите признаки с коэффициентами, близкими к нулю (< 1e-3). Как вы думаете, почему модель исключила их из модели (вновь посмотрите на графики)? Верно ли, что они никак не влияют на спрос на велосипеды?
    • Ответ: windspeed(ms) - этот признак влияет, просто его лассо обрезало регуляризацией в силу мультиколлинеарности, что уже обсуждалось.

In [19]:
zip(df.columns[0:(df.shape[1]-1)], clf.coef_)


Out[19]:
[('season', 532.01898284135359),
 ('yr', 1015.0602226430595),
 ('mnth', -100.03952614356648),
 ('holiday', -83.293959875298938),
 ('weekday', 132.50446549095784),
 ('workingday', 51.557085614073905),
 ('weathersit', -330.55985673998129),
 ('temp', 370.67985503003604),
 ('atemp', 581.39693106549578),
 ('hum', -140.00740550068869),
 ('windspeed(mph)', -191.77140847135024),
 ('windspeed(ms)', -2.685510542868954e-08)]

Заключение

Итак, мы посмотрели, как можно следить за адекватностью линейной модели, как отбирать признаки и как грамотно, по возможности не настраиваясь на какую-то конкретную порцию данных, подбирать коэффициент регуляризации.

Стоит отметить, что с помощью кросс-валидации удобно подбирать лишь небольшое число параметров (1, 2, максимум 3), потому что для каждой допустимой их комбинации нам приходится несколько раз обучать модель, а это времязатратный процесс, особенно если нужно обучаться на больших объемах данных.